Exotic vortex lattices in two-species Bose-Einstein condensates 
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We numerically investigate vortex lattices in rotating two-component Bose-Einstein condensates 
in which the two components have unequal atomic masses and interact attractively with each other. 
For sufficiently strong attraction, the system is found to exhibit exotic ground-state structures in a 
harmonic trap, such as lattices having a square geometry or consisting of two-quantum vortices. The 
obtained states satisfy the Feynman relation, and they can be realized with current experimental 
techniques. 
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I. INTRODUCTION 

One of the characteristic properties of superfluids is 
that they respond to rotation by forming quantized vor- 
tices P, [H • Perhaps the most convenient environment to 
controllably study the physics of vortices is provided by 
dilute Bose-Einstein condensates (BECs) of alkali-metal 
atoms [H-IH . The experiments with dilute BECs have ver- 
ified Abrikosov's original prediction [6j that when a con- 
densate without spin degrees of freedom is set into rapid 
rotation, a triangular lattice of single-quantum vortices 
is created. The areal density of vortices n v is given by 
the so-called Feynman relation [2[ 



(1) 



where m is the mass of the constituent boson and f2 is 
the angular frequency at which the superfluid is being 
rotated. Equation (jTJ) is obtained by assuming that, on 
average, the superfluid rotates like a rigid body — even 
though the superfluid flow is inherently irrotational and 
has a nonvanishing curl only at the singular vortex cores. 

Vortex physics becomes much more diverse in BECs 
that consist of more than one component Q. Even for 
the simplest multicomponent case, the two-component 
BEC, a rich variety of exotic ground-state vortex struc- 
tures have been found in theoretical studies: coreless 
vortices that can be interpreted as skyrmions or meron 
pairs using a pseudospin representation [a, B , serpen- 
tine vortex sheets , giant skyrmions [U [l2| , and in- 
terlacing square vortex lattices [12-fl8|. Of these, the 
singly quantized coreless vortices [3J and the square lat- 
tices [19| have already been observed experimentally. To 
date, two-component BECs have been realized in various 
setups: a single isotope in two different hyperfine spin 
states yJT^fij] , two different isotopes of the same alkali 
metal |24l.l25|. or two distinct elements [2ril - [3lj | . 

All the aforementioned vortex structures pertain to 
systems where the two components interact repulsively 
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with each other. Indeed, the tendency to form structures 
where one component fills the vortex core of the other 
component is understood in terms of the minimization 
of the repulsive density-density interaction energy. Also 
the transition from triangular to square vortex lattices is 
driven by the repulsive intercomponent interaction [l~2l — 

EH- 

However, vortex lattices in rotating two-component 
BECs with attractive intercomponent interaction have 
not been studied as thoroughly (3^, HI] . When all par- 
ticles in the system have equal mass and experience the 
same external rotation, the two components are expected 
to form identical vortex lattices and behave in most re- 
spects like a single-component BEC. If the atomic masses 
of the two components are unequal, the situation be- 
comes more interesting: On one hand, Eq. (p} implies 
that the vortex densities in the two components should 
differ from one another. The intercomponent attraction, 
on the other hand, results in an effective attraction be- 
tween vortices in different components, favoring config- 
urations in which the two vortex lattices have maximal 
overlap. Clearly, differently spaced triangular lattices can 
overlap without frustration only for a very limited range 
of atomic mass ratios. Thus, we pose the following ques- 
tions: What is the combined effect of the mass difference 
and the intercomponent attraction on the arrangement 
of vortices in the rotating two-component system? Is the 
Feynman relation still obeyed? 

In this article, we aim to answer the above questions. 
We demonstrate numerically that the two-component 
BEC with unequal atomic masses and attractive inter- 
component interaction supports exotic ground-state vor- 
tex configurations in a rotating harmonic trap. Depend- 
ing on both the ratio of atomic masses and the strength 
of the intercomponent attraction, the ground states may 
contain overlapping square lattices or even lattices con- 
sisting of doubly quantized vortices. We find no signif- 
icant deviations from the Feynman relation. Since two- 
atomic-species BECs have already been experimentally 
realized 126143 1| . also with tunable intercomponent inter- 
actions [29( , the exotic lattices should be observable with 
the current state of the art. The two-component sys- 
tem with unequal masses was first studied by Barnett et 
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al. [32|, l33|, but they considered only mass ratios near 
unity and thus did not detect any nontriangular vortex 
lattices. 

The remainder of this article is organized as follows. In 
Sec. HH we present the zero-temperature mean-field the- 
ory of the two-component BEC and derive the Thomas- 
Fermi approximation that provides useful information 
about the ground-state structure. We also generalize the 
Feynman relation and discuss the expected lattice geome- 
tries for selected values of the atomic mass ratio. Sec- 
tion IIIII presents the numerically obtained ground-state 
vortex lattices. In Sec. HVI we summarize and discuss the 
main results of the work. 



II. THEORY 

We restrict our studies to the zero-temperature mean- 
field regime, where the two-component BEC is described 
in terms of two complex- valued order-parameter fields 
and ^2- The system is set into rotation about the z axis 
with angular frequency f2 which is assumed to be the 
same for both components. For simplicity, we only con- 
sider situations in which the z dependence of and \&2 
can be factored out and all the relevant dynamics occur 
in the x and y directions. This reduction is accurate, e.g., 
in the regime of strong axial confinement, where the har- 
monic trapping frequencies in the z direction are much 
greater than the trapping frequencies u>\ and uji in the 
radial direction. The two-dimensional Gross-Pitaevskii 
(GP) energy functional in the rotating reference frame is 
given by E = J e(r)d 2 r, where the energy density is 



N2, yields the coupled GP equations 



ft 2 V 2 
2m i 



+ l m jU 2 r 2 + 9jnj 
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%(r)=0, 



(3) 



where j £ {1,2}. Equations ([3]) form the backbone of 
our numerical analysis in Sec. IIIII 



A. Thomas— Fermi approximation 

Let us consider the Thomas-Fermi (TF) limit corre- 
sponding to large particle numbers [3 a . |36j, where the 
spatial variation of the condensate density is negligi- 
ble compared to the other terms in the energy func- 
tional. Thus, we approximate — ifiV^j = HCVSj^j 



ihexp (iSj) V\^j 



jVj'fyj, where Sj — arg(\Pj) and 



Vj = hVSj/rrij denotes the superfluid velocity of compo- 
nent j. Consequently, Eq. ^ simplifies to 

1 2 

^tfO) = - ^2 [ m i ( V J _ Vl ' b ) 2 n i + mjtfr 2 nj 



+9j n j I + 3l2™l™2, 



(4) 



where u>j := uj 2 — Q 2 and v r b = Qe z x r corresponds to 

classical rigid-body rotation. For large BECs with many 
vortices, the superfluid flow closely mimics the rigid-body 
behavior, that is, Vj ~ v r b outside the vortex cores, and 
the first term vanishes to leading order. Hence, we end 
up with 



2 1 2 

e(r) - J2 ( — |V*i| 2 + - mj uj 2 r 2 nj + -g 3 n 2 e '( r ) = 3 £ \ m ^V n 3 + 9jn 2 ] + g^mm, 



(5) 
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Here, m\ and mi are the masses of the bosons in the two 



components, r 



y , and L z 



-ih [e z ■ (r x V)] 



is the z component of the angular momentum operator 
in the position representation. The order parameters are 
normalized such that J nj(r)d 2 r — Nj, where rij = \^j\ 2 
is the areal particle density and Nj the total particle 
number of component j € {1,2}. In Eq. @, gi and g%, 
which measure the intracomponent interaction strengths, 
are assumed to be positive throughout the paper, and 1712 
is the intercomponent coupling constant for which we 
consider only negative values. The parameter gj (1712) 
is proportional to the corresponding s-wave scattering 
length dj (012), with the exact dependence being deter- 
mined by the trap along the z direction [34]]. 

Minimizing the GP energy functional with respect to 
variations of 'Si and ^2, and introducing the chemical 
potentials /ii and /12 to fix the particle numbers Ni and 



which has exactly the same form as the TF energy func- 
tional for a nonrotated two-component system, but the 
confining potentials are altered to rrij (lu 2 — il 2 ) r 2 /2. 

Variation of s' with respect to n\ and 712 with fixed par- 
ticle numbers J rijd 2 r — Nj gives two coupled equations 
that can be solved for n\ and ni. When g\gi — g\ 2 > 0, 
the components are miscible, and we find 



nj(r)=nj !0 (l-r 2 /R 2 j), 



(6) 



where the peak densities {nj t o} and TF radii {Rj} are 
given by 



Nj gz-jmjU] 2 - g 12 ms-joj 2 ^j 
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5i 52 - 9l2 



n 9a-jmjQ? - 9\im?,-ju\_j 



(7) 
(8) 



The healing lengths, which can be used to estimate the 



vortex core radii, are defined as £j = hj {2mjji'jj 
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where p'j — gjUj.o + 512^3-7,0 denotes the chemical po- 
tential in the TF approximation. 

In order to simplify the numerical analysis, it is conve- 
nient reduce the number of free parameters in the prob- 
lem. Thus, we will fix the parameters pertaining to com- 
ponent 2, i.e., w r .2, 52, and N%, in a way that is most 
suitable for studying how the vortices reorganize when 
the interaction strength 512 is varied. To this end, we 
require the TF profiles to be identical, n^o = «2,o = : n a 
and R\ = R2 =: R, as well as equating the healing 
lengths, & = £ 2 =: £• According to Eqs. © and ©, 
these conditions are met when 



n 2 + p 2 (u 



2 

92 = P9i + {p 
N 2 =N 1 , 



512, 



(9) 



where p = mi/1712 denotes the atomic mass ratio. Con- 
sequently, the TF expressions simplify to 



TOlATi lo 2 - tt 2 



7T 51+512 
r4 = 4N\ g x +912 



Trh 4 



n 2 



(10) 
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We note that the interaction strengths enter Eqs. ([TO]) 
only through the combination g\ + 512 which may be in- 
terpreted as measuring the total density-density coupling 
strength. Therefore, it will be convenient to parametrize 
the system in terms of gi+gi2 and the ratio a := —512/51, 
where < a < 1. The case a > 1 corresponds to the col- 
lapse of the BEC under attractive forces. 



B. Feynman relation and vortex lattices 

For two-component BECs, the Feynman relation ([1]) 
naturally generalizes to 



rrijQ 



(11) 



It states that for a two-component system with unequal 
masses, the vortex densities are also unequal, satisfying 
7i v .i = jdn v ,2- Without loss of generality, we will assume 
that p > 1. It should be noted that Eq. (1TT1) differs signif- 
icantly from the treatment of Barnett et al. [32|, [33[ who 
equate the two vortex densities and allow for component- 
dependent rotation frequencies. 

How does the difference in the vortex densities affect 
the geometry of the underlying lattice? The attractive in- 
tercomponent interaction favors states in which vortices 
in different components lie on top of each other, since it 
acts to maximize the overlap J niri2d 2 r between the two 
components. However, two regular triangular lattices of 
single-quantum vortices can perfectly overlap only when 
the mass ratio can be expressed as p — k 2 + I 2 + kl, 



where k, I S N. Therefore, the two-component BEC may 
support unconventional lattice geometries when p is far 
from such a value and the intercomponent attraction is 
sufficiently strong. 

The overlapping vortex lattices are illustrated in Fig. [I] 
For simplicity, we consider here only single-quantum vor- 
tices |37j and integer values of p. In real experiments, 
p is rarely strictly integer-valued but often sufficiently 
close to an integer such that the lattices in Fig. [1] pro- 
vide good approximations to the actual ones. For p = 1 
and 512 < 0, the ground state consists of two superposed 
triangular lattices [Fig. 1(a)]. In the experimentally rel- 
evant case p = 2 (for 41 K- 87 Rb BEC (M^, p « 2.1), 
a natural way to obtain n Vj i/n v _2 = 2 is to arrange the 
vortices into two square lattices tilted by 45 degrees rel- 
ative to one another [Fig.[IJb)]. For p = 3, which applies 
to 41 K- 133 Cs (p w 3.2) and 7 Li- 23 Na (p « 3.3), triangu- 
lar lattices are expected [Fig. Hfc)]. This is also the case 
for p — 4 (in 23 Li- 87 Rb mixtures, p ~ 3.8), as shown 
in Fig. Hfd). It is also possible to superpose two square 
lattices with n Vi i/n Vi 2 — 4 [Fig. [Tfe)], but such a state 
is expected to have a higher energy than the triangular 
arrangement due to the higher intracomponent energies. 
In the next section, we present numerical solutions of the 
two-component GP equations which exhibit the lattice 
geometries depicted in Fig. Q] 



III. NUMERICAL RESULTS 

In order to demonstrate that the intercomponent inter- 
action can induce a change in the ground-state vortex ge- 
ometry, we have numerically solved the coupled GP equa- 
tions of a rotating two-component BEC with mi 7^ 1712. 
We consider the mass ratios p — m 1 /m 2 £ {2,3,4}. In 
the numerical simulations, Eqs. ([3]) are discretized on 
a uniform grid and solved using a relaxation method. 
We measure length in units of the radial harmonic os- 
cillator length a r — y/ h/miuii and energy in units of 
fuui and normalize the dimensionless order parameters 
to unity. The interactions are parametrized by the di- 
mensionless quantities 5tot := (51 + 512) ^l/fi^ics 2 . and 
a := —512/51- The values we choose for the total inter- 
action parameter g tot are large enough such that the TF 
expressions, Eqs. ([6J and (|10p . accurately approximate 
the numerically obtained states, thereby justifying our 
use of Eqs. ©. To emphasize the role of the intercom- 
ponent attraction, we consider a range of values for a in 
the interval < a < 1. In experiments, the ratio 512/51 
may be controlled with intercomponent Feshbach reso- 
nances [38j], which have been demonstrated for various 
alkali- metal mixtures [2^, [39l - l4ll ] . 

Figure [2ja) presents the ground-state density profiles 
n\ and 712 for the mass ratio p — 2, rotation frequency 
O = 0.97 cji, total interaction parameter g to t = 705, and 
different values of the intercomponent attraction strength 
a . At a = 0, the components are noninteracting, and the 
ground state contains triangular lattices in both com- 
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FIG. 1. (Color online) Vortex lattice geometries of the attrac- 
tively interacting two-component BEC for different values of 
the mass ratio p = mi/m2; (a) p = 1; (b) p = 2; (c) p = 3; 
(d)-(e) p = 4. Hollow and solid circles represent vortices of 
components 1 and 2, respectively. The unit cell of the com- 
bined lattice is shown with dashed lines. 



ponents. However, as the intercomponent attraction in- 
creases, the ground-state profiles change gradually into 
square lattices as shown for a — 0.6 [see also Fig. QJb)]. 
A further increase in a eventually results in the trian- 
gular geometry presented in Fig. [lja), but with vortex 
pairs occupying the lattice points in component 1 instead 
of solitary vortices. The vortex pairs in overlap with 
single vortices in ^ that have elongated cores. Finally, 
at a = 0.97, the vortex pairs have merged to form a tri- 
angular lattice of two-quantum vortices in ^i, while ^>2 
contains an essentially identical arrangement of single- 
quantum vortices. Although we present results only for 
the values of f2 and gtot mentioned above, qualitatively 
similar behavior is observed also for other values as long 
as the states contain a sufficient number of vortices. We 
point out that the square vortex lattice shown here for 
a = 0.6 also appears in a system of two homogeneous 
superfluids with optically induced current-current inter- 
actions and p = 2 [42j |. 

To assess the validity of the Feynman relation, 
Eq. (fTTj) . it is useful to investigate the velocity fields 
vj = hVSj/rrij. We define the average azimuthal ve- 



0.00 



0.60 



0.75 



0.97 



(b) 


2.5r 




2.0- 









1.5* 


'5 




l£ 


1.0- 


IN 

IS 


0.5- 


IS 






FIG. 2. (Color online) Ground states of a two-component 
BEC with the mass ratio mi/m.2 = 2, rotation frequency 
Q, — 0.97 tJi, and interaction strength g to t = 705. (a) Particle 
densities n\ and 722 at different values of the intercomponent 
attraction strength a = — 512/(71 ■ (b) Average azimuthal ve- 
locities vi and V2 as functions of r [Eq. (|12[) ] for the states 
shown in (a). For clarity, the curves have been shifted ver- 
tically; the dotted line TT r b = £lr passes through the origin 
in all c ases. Th e Thomas-Fermi radius R — 11.1 a r , where 
a r = \Jhjvri\Lki\. 
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locities as 



(a) 



Vj{r) 



Jo* e <f> ■ v j( r ' ^K( r ' <t>) d( t> 



Jo n 3 



V, cj))d(j) 



(12) 



where 4> is the polar angle. The average is weighted with 
the particle density to smooth out the velocity diver- 
gences at the singular vortex cores. For classical rigid- 
body rotation, the velocity field is v r b = £le z x r, and 
hence Eq. (fl~2|) gives v T b(r) = Vtr. If Eq. (fTTj) holds and 
the supcrfluid flow mimics the rigid-body behavior, v\{r) 
and V2(r) should approximately fall on top of the line fir. 
In Fig. [2fb), the velocities are plotted as functions of r 
for the states shown in Fig. [2ja). We observe that v\ 
and T>2 closely follow the line fir, in accordance with the 
Feynman relation, irrespective of the value of a. When 
cr is sufficiently large, v\ and T>2 lie accurately on top of 
each other, reflecting the fact that the vortex cores in the 
two components overlap. 

For the mass ratio p = 3, two triangular lattices can 
efficiently overlap [Fig. QJc)], and consequently we have 
not observed any nontriangular geometries for large lat- 
tices, despite varying the parameters ft, gtot, and a ex- 
tensively. Representative ground states for p = 3 are 
depicted in Fig. [3ja) . The vortex lattices retain their tri- 
angular shape with increasing a. However, those vortices 
in that are paired with vortices in ^2 increase in core 
size as a increases, whereas unpaired vortices become 
smaller. A comparison of the various vortex core sizes 
with the healing length £ is presented in Fig. EJc) , where 
the core size at a = is observed to lie between the core 
sizes of the paired and unpaired vortices at a — 0.975. 
Moreover, the unpaired vortices of component 1 are found 
to induce density craters in component 2 at sufficiently 
large values of the intercomponent attraction as shown in 
Fig.^a) for a = 0.975. The azimuthal velocities plotted 
in Fig. [3jb) again show no significant deviation from the 
Feynman relation. 

For p = 4, the ground state is found to consist of 
two triangular lattices as expected [Fig. [TJd)]. The nu- 
merical results are exemplified in Fig. Eta), where the 
ground-state particle densities are presented for p = 4, 
£1 = 0.97 wi, <7tot = 4000, and different values of a. 
Paired vortices are again found to have larger cores than 
unpaired ones, and craters develop in ri2 at the positions 
of unpaired vortices in at large values of a. Plots of 
the average velocities Vj [Fig. HJb)] faithfully reproduce 
the rigid-body behavior in accordance with the Feynman 
relation. 

Although the triangular lattice is always the ground 
state for p = 4, its energy at large a turns out to be 
nearly degenerate with a stationary square-lattice state. 
An example of such a state is shown in Fig. [4jc) for 
a = 0.85, and it corresponds to the lattice geometry 
illustrated in Fig. [TJe). The relative energy difference 
between this state and the ground state presented in 
Fig. 0|a) is AE/E k 4x 10~ 5 . Various metastable 
structures may thus appear in this parameter region. To 
study whether the state in Fig. @|c) is robust enough to 
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FIG. 3. (Color online) Ground states of a two-component 
BEC with mi/m 2 = 3, CI = 0.95 ui, and g tot = 2000. 

(a) Densities m and 712 at different values of a — —312/51. 

(b) Average azimuthal velocities vi and U2 [Eq. (I12|l ] for the 
states presented in (a). The curves have been shifted in the 
vertical direction for clarity, and R — 12.7 a r . (c) Magnified 
view of ni for —512/31 = and 0.975 [boxed areas in (a)] 
illustrating the variation of the vortex core size. The circles 
have the radius £ = 0.252 a r . 



be experimentally realizable, one could solve its elemen- 
tary excitation spectrum. However, this investigation is 
left for future work. 
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FIG. 4. (Color online) Vortex lattices for roi/m2 = 4, f2 = 
0.97 cji, and gtot = 4000. (a) Ground-state densities ni and 
7i2 at different values of a. (b) Average velocities v\ and 1J2 
for the states presented in (a). The curves have been shifted 
vertically, and R = 17.1 a r . (c) Stationary state at a = 0.85 
having only slightly higher energy than the ground state. 



vortex configurations, such as square and two-quantum- 
vortex lattices, when the intercomponent attraction is 
sufficiently strong and the ratio of atomic masses is suit- 
able. Importantly, the nontriangular geometry or the 
multiquantum nature of the lattices was not induced by 
external fields [43[ but emerged as an inherent property 
of the system. The pairing of vortices between the two 
components was also observed to affect the vortex core 
size, paired vortices being larger than unpaired ones. 

We also investigated the validity of the Feynman rela- 
tion in the two-component system and found no signifi- 
cant deviations from it. This should be contrasted with 
the findings of Barnett et al. who studied the attrac- 
tively interacting two-component BEC for mass ratios 
only slightly above unity, 1 < p < 1.6, and concluded 
that interlocking of the vortices can violate the Feynman 
relation. However, this anomaly occurred only within 
a finite distance from the trap center, whereas in the 
states considered in this work the interlocking behavior 
extended through the whole system. Therefore, interest- 
ing crossover behavior in the regime 1 < p < 2 may occur; 
this would be relevant to 87 Rb- 133 Cs mixtures [3lT l4lj|. 
for which p = 1.5. 

Perhaps the most promising system for observing the 
nontriangular vortex lattices is the two-species BEC of 
41 K and 87 Rb, which has already been realized with tun- 
able intercomponent interactions 29] . This combination 
yields the mass ratio p — 2.1, and we have verified nu- 
merically that such a system exhibits ground-state vortex 
structures similar to those found for p = 2 [Fig. [5] . In 
particular, the two-quantum-vortex lattice expected for 
strong intercomponent attraction would be a rare exam- 
ple of a ground state containing multiply quantized vor- 
tices with self-supporting, genuinely empty cores. Since 
the energy of a vortex increases quadratically with its 
quantum number, multiply quantized vortices do not 
usually appear in the ground state. 

There are many ways to extend this work in the future. 
For example, it would be interesting to study the system 
in the limit of very fast rotation by generalizing the an- 
alytical calculations of Mueller and Ho [l3| to dissimilar 
condensates (cf. also Refs. [H, EH). In the simulations, 
one could also scan different parameter values more ex- 
tensively to obtain vortex phase diagrams in the (cr, fi) 
plane [H, [l4| . This could be done efficiently by formu- 
lating the ener geti cs in terms of the positions of vortices 
as in Refs. \MM^- 



ACKNOWLEDGMENTS 



IV. DISCUSSION 

In summary, we have studied vortex lattices in rotat- 
ing two-component BECs with an atomic mass difference 
and attractive interaction between the two components. 
We found that such systems support exotic ground-state 



We gratefully acknowledge financial support from the 
Academy of Finland, the Emil Aaltonen Foundation, the 
Finnish Cultural Foundation, the KAUTE Foundation, 
and the Vaisala Foundation. We thank E. Ruokokoski 
for insightful discussions and V. Pictila and S. M. M. Vir- 
tanen for their help during the early stages of the work. 



7 



L. Onsager, Nuovo Cimento 6, Suppl. 2, 249 (1949). 
R. P. Feynman, in Progress in Low Temperature Physics, 
edited by C. J. Gorter (North-Holland Publishing Com- 
pany, Amsterdam, 1955), Chap. 2 

M. R. Matthews, B. P. Anderson, P. C. Hal jan, D. S. 
Hall, C. E. Wieman, and E. A. Cornell, [Phys. Rev. Lett.| 
83, 2498 (1999). 

K. W. Madis on, F. Ch evy, W. Wohlleben, and J. Dal- 
ibard, [Phys! Rev. Lett.| 84, 806 (2000). 
J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ket- 
terle, Science 292, 476 (2001). 

A. A. Abrikosov, Zh. Eksp. Teor. Fiz. 32, 1442 (1957) 
[Sov. Phys. JETP 5, 1174 (1957)] 

K. Kasamatsu, M. Tsubota, and M. Ueda, Int. J. Mod. 
Phys. B 19, 1835 (2005). 

K. Kasamatsu, M. Tsubota, and M. Ueda, 

|Phys. Rev. LetE] 93, 250406 (2004). 

K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. A 
71, 043611 (2005). 

K. Kasamatsu and M. Tsubota, [Fhys. Rev. A| 79, 023606 
(2009). 

S.-J. Yang, Q.-S. Wu, S.-N. Zhang, and S. Feng, 
|Phys. Rev. A] 77, 033621 (2008). 

P. Mason and A. Aftalion, Phys. Rev. A 84, 033611 
(2011). 

E. J. Mueller and T.-L. Ho, Phys. Rev. Lett. 88, 180403 
(2002). 

K. Kasamatsu, M. Tsubota, and M. Ueda, 

|Phys. Rev. Lett~| 91, 150406 (2003). 

M. Keceli and M. 6. Oktel, Phys. Rev. A 73, 023611 

(2006). 

M. P. Mink, C. M. Smith, and R. A. Duine, Phys. Rev. 
A 79, 013605 (2009). 

J. A. M. Huhtam&ki, T. P. Simula, M. Kobayashi, and 
K. Machida, Phys. Rev. A 80, 051601 (2009). 
E. K. Dahl, E. Babaev, and A. Sudb0, Phys. Rev. B 78, 
144510 (2008). 

V. Schweikhard, I. Coddington, P. Engels, S. Tung, and 
E. A. Cornell, [Phys. Rev. Lett"] 93, 210403 (2004). 

C. J. Myatt, E. A. Burt, R. W. Ch rist, E. A. Cornell, 
and C. E. Wieman, phys. Rev. Lettl. 78, 586 (1997). 

D. S. Hall, M. R. Matthews, J. R. En sher, C. E. Wieman, 
and E. A. Cornell, [Phys. Rev. LetE] 81, 1539 (1998). 

G. Delannoy, S. G. Murdoch, V. Boyer, V. Josse, 
P. Bouyer, and A. Aspect, |Phys. Rev. "A] 63, 051602 
(2001). 

R. P. Anderson, C. Ticknor, A. I. Sidorov, and B. V. 

Hall, |Phys. Rev. A| |80, 023603 (2009). 

S. B. Papp, J. M. Pino, and C. E. Wieman, 

|Phys. Rev. Lettl lOl, 040402 (2008). 

S. Sugawa, R. Yamazaki, S. Taie, and Y. Takahashi, 



[26 

[27 
[28 
[29 

[30 

[31 

[32 
[33 
[34 



[35 
[36 

[37 



[38 
[39 

[40 
[41 

[42 
[43 
[44 



Phys. Rev. A] 84, 011610 (2011). 

G. Ferrari, M. Inguscio, W. Jastrzebski, G. Modugno, 
G. Roati, and A. Simoni, |Phys. Rev. Lett.| 89, 053202 
(2002). 

G. Modugno, M. Modugno, F. Riboli, G. Roati, and 

M. Inguscio, Phys. Rev. Lett. 89, 190404 (2002). 

J. Catani, L. De Sarlo, G. Barontini, F. Minardi, and 

M. Inguscio, Phys. Rev. A 77, 011603 (2008). 

G. Thalhammer, G. Barontini, L. De Sarlo, J. Catani, 

F. Minardi, and M. Inguscio, |Phys. Rev. Lett.| 100, 

210402 (2008). 

K. Aikawa, D. Akamatsu, J. Kobayashi, M. Ueda, 
T. Kishimoto, and S. Inouye, New J. Phys. 11, 055035 
(2009). 

A. Lercher, T. Takekoshi, M. Debatin, B. Schuster, 
R. Rameshan, F. Ferlaino, R. Grimm, and H.-C. Nagerl, 
Eur. Phys. J. D 65, 3 (2011). 

R. Barnett, G. Refael, M. A. Porter, and H. P. Biichler, 
New J. Phys. 10, 043030 (2008). 

R. Barnett, E. Chen, and G. Refael, New J. Phys. 12, 
043004 (2010). 

In the case of strong harmonic confinement in the z 
direction with trapping frequencies {uu Zt j}, the cou- 
pling constants are given by gj 



/8nh~ aj /m,ja Zt j 
where a z , = 



1/2 



and gi2 = 012/ 'm\i 
^/h/mjLUzj and mi2 = mwa^l (mi + 7712) is the reduced 
mass. 

A. L. Fetter, Rev. Mod. Phys. 81, 647 (2009). 

I. Corro, R. G. Scott, and A. M. Martin, Phys. Rev. A 

80, 033609 (2009). 

Vortices with k > 1 quanta are not expected in the 
ground state, because their energy scales as k 2 . Neverthe- 
less, two-quantum vortices are observed for p = 2 when 
the intercomponent attraction is very strong (Sec. IIII| ). 
C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, 
|Rev. Mod. Phys1, 82, 1225 (2010). 

C. A. Stan, M. W. Zwierle in, C. H. Schunc k, S. M. F. 
Raupach, and W. Ketterle, [Phys. Rev. Lett] 93, 143001 
(2004). 

S. Inouye, J. Goldwin, M. L. Olsen, C. Ticknor, J. L. 
Bohn, and D. S. Jin, |Phys. Rev. Lett.| 93, 183201 (2004). 
K. Pilch, A. D. Lange, A. Prantner, G. Kerner, F. Fer- 
laino, H.-C. Nagerl, and R. Grimm, Phys. Rev. A 79, 
042718 (2009). 

E. K. Dahl, E. Babaev, and A. Sudb0, Phys. Rev. Lett. 
101, 255301 (2008). 

S. Tung, V. Schweikhard, and E. A. Cornell, 
|Phys. Rev. LetE] 97, 240402 (2006). 
A. Aftalion, P. Mason, and J. Wei, Phys. Rev. A 85, 
033614 (2012). 



